function [ newFovRow newFovCol ] = findMpMapCenter( ratioIm, imageFeatures )
%findMpMapCenter finds the macular pigment center
%   Detailed explanation goes here

% Define region close to the original center to find minimum pixels in.
% This is done to avoid the selection of minimum pixels along the macular
% boundary (see the macular pigment maps and you will notice the black dots
% outlining the circular region).
BOXSIZE = round(imageFeatures.macRadius/1.5);
rowRange = imageFeatures.fovRow - BOXSIZE: imageFeatures.fovRow + BOXSIZE;
colRange = imageFeatures.fovCol - BOXSIZE: imageFeatures.fovCol + BOXSIZE;
validPixelValues = ratioIm(rowRange, colRange);
% Determine minimum pixel value in the region of interest. Take the
% absolute value to avoid negative values, this value will be used in
% thresholding.
minPixValue = min(min(abs(validPixelValues)));
for percentValues = [0.15 0.13 0.10 0.08 0.05 0.03 0.01];
    threshold = minPixValue + percentValues*minPixValue;
    if threshold < 1
        break;
    end
end

% Binary the image according to a threshold. The white pixels will all be
% potential choices for the new fovea center
binaryRatioIm = ~im2bw(ratioIm, threshold);

% Occasionally dark pixels are found at the rim of the macular boundary due
% to the image processing of dividing and normalizing (see the macular
% pigment maps and you will notice the black dots outlining the circular
% region). We need to set these pixel values outside the valid region equal
% to 0 since valid values in the binary image are equal to 1.

maskIm = zeros(size(ratioIm));
maskIm(rowRange, colRange) = 1;
validBinaryRatioIm = immultiply(binaryRatioIm, maskIm);
potentialCenterPixels = find(validBinaryRatioIm == 1);
[potentialRows potentialCols] = ind2sub(size(validBinaryRatioIm), potentialCenterPixels);
newFovRow = round(mean(potentialRows));
newFovCol = round(mean(potentialCols));

end

